!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine bnd_grpr (ro, mx, my, mz, bx, by, bz, en, gx, gy, gz, x, y, z, gm, nx, ny, nz, ng, dir, sides)
!==================================================================================================================================|
!  sides = 'l' for lower boundary
!  sides = 'h' for higher boundary
!  sides = 'b' for both
!----------------------------------------------------------------------------------------------------------------------------------|

   implicit none

   character(1), intent(in) ::  dir, sides
   integer(4), intent(in) :: ng, nx, ny, nz
   real(8), intent(in) :: gm
   real(8) :: x(nx), y(ny), z(nz)
   real(8), dimension(nx, ny, nz), intent(inout) :: ro, mx, my, mz, bx, by, bz, en, gx, gy, gz

   integer(4) :: i, j, k, ngp1, ngp2, nxm1, nxm0, nym1, nym0, nzm1, nzm0, bnd_1, bnd_2
   real(8) :: pr, v2, b2, gmm1
   
!----------------------------------------------------------------------------------------------------------------------------------|
   gmm1 = gm - 1.0d0
   ngp1 = 1 + ng
   ngp2 = 2 + ng

   if (dir .eq. 'x') then

      if (sides .eq. 'l' .or. sides .eq. 'b') then
         do k = 1, nz
         do j = 1, ny
         do i = 1, ng
            bnd_2 = ngp1 - i
            bnd_1 = ngp2 - i
            v2 = mx(bnd_1, j, k) * mx(bnd_1, j, k) + my(bnd_1, j, k) * my(bnd_1, j, k) + mz(bnd_1, j, k) * mz(bnd_1, j, k)
            b2 = bx(bnd_1, j, k) * bx(bnd_1, j, k) + by(bnd_1, j, k) * by(bnd_1, j, k) + bz(bnd_1, j, k) * bz(bnd_1, j, k)
	         pr = (en(bnd_1, j, k) - v2 / ro(bnd_1, j, k) / 2.0d0 - b2 / 2.0d0) * gmm1
            pr = pr + (ro(bnd_2, j, k) * gx(bnd_2, j, k) + ro(bnd_1, j, k) * gx(bnd_1, j, k)) / 2.0d0 * (x(bnd_2) - x(bnd_1))
            v2 = mx(bnd_2, j, k) * mx(bnd_2, j, k) + my(bnd_2, j, k) * my(bnd_2, j, k) + mz(bnd_2, j, k) * mz(bnd_2, j, k)
            b2 = bx(bnd_2, j, k) * bx(bnd_2, j, k) + by(bnd_2, j, k) * by(bnd_2, j, k) + bz(bnd_2, j, k) * bz(bnd_2, j, k)
            en(bnd_2, j, k) = v2 / ro(bnd_2, j, k) / 2.0d0 + b2 / 2.0d0 + pr / gmm1
         enddo
         enddo
         enddo
      endif

      nxm0 = nx - ng
      nxm1 = nx - ng - 1

      if (sides .eq. 'h' .or. sides .eq. 'b') then
         do k = 1, nz
         do j = 1, ny
         do i = 1, ng
            bnd_2 = nxm0 + i
            bnd_1 = nxm1 + i
            v2 = mx(bnd_1, j, k) * mx(bnd_1, j, k) + my(bnd_1, j, k) * my(bnd_1, j, k) + mz(bnd_1, j, k) * mz(bnd_1, j, k)
            b2 = bx(bnd_1, j, k) * bx(bnd_1, j, k) + by(bnd_1, j, k) * by(bnd_1, j, k) + bz(bnd_1, j, k) * bz(bnd_1, j, k)
	         pr = (en(bnd_1, j, k) - v2 / ro(bnd_1, j, k) / 2.0d0 - b2 / 2.0d0) * gmm1
            pr = pr + (ro(bnd_2, j, k) * gx(bnd_2, j, k) + ro(bnd_1, j, k) * gx(bnd_1, j, k)) / 2.0d0 * (x(bnd_2) - x(bnd_1))
            v2 = mx(bnd_2, j, k) * mx(bnd_2, j, k) + my(bnd_2, j, k) * my(bnd_2, j, k) + mz(bnd_2, j, k) * mz(bnd_2, j, k)
            b2 = bx(bnd_2, j, k) * bx(bnd_2, j, k) + by(bnd_2, j, k) * by(bnd_2, j, k) + bz(bnd_2, j, k) * bz(bnd_2, j, k)
            en(bnd_2, j, k) = v2 / ro(bnd_2, j, k) / 2.0d0 + b2 / 2.0d0 + pr / gmm1
         enddo
         enddo
         enddo
      endif

   elseif (dir .eq. 'y') then

      if (sides .eq. 'l' .or. sides .eq. 'b') then
         do k = 1, nz
         do j = 1, ng
            bnd_2 = ngp1 - j
            bnd_1 = ngp2 - j
            do i = 1, nx
               v2 = mx(i, bnd_1, k) * mx(i, bnd_1, k) + my(i, bnd_1, k) * my(i, bnd_1, k) + mz(i, bnd_1, k) * mz(i, bnd_1, k)
               b2 = bx(i, bnd_1, k) * bx(i, bnd_1, k) + by(i, bnd_1, k) * by(i, bnd_1, k) + bz(i, bnd_1, k) * bz(i, bnd_1, k)
	            pr = (en(i, bnd_1, k) - v2 / ro(i, bnd_1, k) / 2.0d0 - b2 / 2.0d0) * gmm1
               pr = pr + (ro(i, bnd_2, k) * gy(i, bnd_2, k) + ro(i, bnd_1, k) * gy(i, bnd_1, k)) / 2.0d0 * (y(bnd_2) - y(bnd_1))
               v2 = mx(i, bnd_2, k) * mx(i, bnd_2, k) + my(i, bnd_2, k) * my(i, bnd_2, k) + mz(i, bnd_2, k) * mz(i, bnd_2, k)
               b2 = bx(i, bnd_2, k) * bx(i, bnd_2, k) + by(i, bnd_2, k) * by(i, bnd_2, k) + bz(i, bnd_2, k) * bz(i, bnd_2, k)
               en(i, bnd_2, k) = v2 / ro(i, bnd_2, k) / 2.0d0 + b2 / 2.0d0 + pr / gmm1
            enddo
         enddo
         enddo
      endif

      nym0 = ny - ng
      nym1 = ny - ng - 1

      if (sides .eq. 'h' .or. sides .eq. 'b') then
         do k = 1, nz
         do j = 1, ng
            bnd_2 = nym0 + j
            bnd_1 = nym1 + j
            do i = 1, nx
               v2 = mx(i, bnd_1, k) * mx(i, bnd_1, k) + my(i, bnd_1, k) * my(i, bnd_1, k) + mz(i, bnd_1, k) * mz(i, bnd_1, k)
               b2 = bx(i, bnd_1, k) * bx(i, bnd_1, k) + by(i, bnd_1, k) * by(i, bnd_1, k) + bz(i, bnd_1, k) * bz(i, bnd_1, k)
	            pr = (en(i, bnd_1, k) - v2 / ro(i, bnd_1, k) / 2.0d0 - b2 / 2.0d0) * gmm1
               pr = pr + (ro(i, bnd_2, k) * gy(i, bnd_2, k) + ro(i, bnd_1, k) * gy(i, bnd_1, k)) / 2.0d0 * (y(bnd_2) - y(bnd_1))
               v2 = mx(i, bnd_2, k) * mx(i, bnd_2, k) + my(i, bnd_2, k) * my(i, bnd_2, k) + mz(i, bnd_2, k) * mz(i, bnd_2, k)
               b2 = bx(i, bnd_2, k) * bx(i, bnd_2, k) + by(i, bnd_2, k) * by(i, bnd_2, k) + bz(i, bnd_2, k) * bz(i, bnd_2, k)
               en(i, bnd_2, k) = v2 / ro(i, bnd_2, k) / 2.0d0 + b2 / 2.0d0 + pr / gmm1
            enddo
         enddo
         enddo
      endif

   elseif (dir .eq. 'z') then

      if (sides .eq. 'l' .or. sides .eq. 'b') then
         do k = 1, ng
            bnd_2 = ngp1 - k
            bnd_1 = ngp2 - k
            do j = 1, ny
            do i = 1, nx
               v2 = mx(i, j, bnd_1) * mx(i, j, bnd_1) + my(i, j, bnd_1) * my(i, j, bnd_1) + mz(i, j, bnd_1) * mz(i, j, bnd_1)
               b2 = bx(i, j, bnd_1) * bx(i, j, bnd_1) + by(i, j, bnd_1) * by(i, j, bnd_1) + bz(i, j, bnd_1) * bz(i, j, bnd_1)
               pr = (en(i, j, bnd_1) - v2 / ro(i, j, bnd_1) / 2.0d0 - b2 / 2.0d0) * gmm1
               pr = pr + (ro(i, j, bnd_2) * gz(i, j, bnd_2) + ro(i, j, bnd_1) * gz(i, j, bnd_1)) / 2.0d0 * (z(bnd_2) - z(bnd_1))
               v2 = mx(i, j, bnd_2) * mx(i, j, bnd_2) + my(i, j, bnd_2) * my(i, j, bnd_2) + mz(i, j, bnd_2) * mz(i, j, bnd_2)
               b2 = bx(i, j, bnd_2) * bx(i, j, bnd_2) + by(i, j, bnd_2) * by(i, j, bnd_2) + bz(i, j, bnd_2) * bz(i, j, bnd_2)
               en(i, j, bnd_2) = v2 / ro(i, j, bnd_2) / 2.0d0 + b2 / 2.0d0 + pr / gmm1
            enddo
            enddo
         enddo
      endif

      nzm0 = nz - ng
      nzm1 = nz - ng - 1

      if (sides .eq. 'h' .or. sides .eq. 'b') then
         do k = 1, ng
            bnd_2 = nzm0 + k
            bnd_1 = nzm1 + k
            do j = 1, ny
            do i = 1, nx
               v2 = mx(i, j, bnd_1) * mx(i, j, bnd_1) + my(i, j, bnd_1) * my(i, j, bnd_1) + mz(i, j, bnd_1) * mz(i, j, bnd_1)
               b2 = bx(i, j, bnd_1) * bx(i, j, bnd_1) + by(i, j, bnd_1) * by(i, j, bnd_1) + bz(i, j, bnd_1) * bz(i, j, bnd_1)
               pr = (en(i, j, bnd_1) - v2 / ro(i, j, bnd_1) / 2.0d0 - b2 / 2.0d0) * gmm1
               pr = pr + (ro(i, j, bnd_2) * gz(i, j, bnd_2) + ro(i, j, bnd_1) * gz(i, j, bnd_1)) / 2.0d0 * (z(bnd_2) - z(bnd_1))
               v2 = mx(i, j, bnd_2) * mx(i, j, bnd_2) + my(i, j, bnd_2) * my(i, j, bnd_2) + mz(i, j, bnd_2) * mz(i, j, bnd_2)
               b2 = bx(i, j, bnd_2) * bx(i, j, bnd_2) + by(i, j, bnd_2) * by(i, j, bnd_2) + bz(i, j, bnd_2) * bz(i, j, bnd_2)
               en(i, j, bnd_2) = v2 / ro(i, j, bnd_2) / 2.0d0 + b2 / 2.0d0 + pr / gmm1
            enddo
            enddo
         enddo
      endif

   else

      write (*, *) 'Direction error from bnd_grpr.'
      stop

   endif

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine bnd_grpr
